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Abstract 

O^ We present a class of high order finite volume schemes for the solution of non- 

|7 j conservative hyperbolic systems that combines the one-step ADER-WENO 

r^ finite volume approach with space-time adaptive mesh refinement (AMR). 

"j^ The resulting algorithm, which is particularly well suited for the treatment 

C of material interfaces in compressible multi-phase fiows, is based on: (i) high 

'— ' order of accuracy in space obtained through WENO reconstruction, (ii) a 

,-H high order one-step time discretization via a local space-time discontinuous 

^ Galerkin predictor method, and (iii) the use of a path conservative scheme for 

(^ handling the non-conservative terms of the equations. The AMR property 

[^ with time accurate local time stepping, which has been treated according 

to a cell-by-cell strategy, strongly relies on the high order one-step time 

discretization, which naturally allows a high order accurate and consistent 

cn computation of the jump terms at interfaces between elements using different 

time steps. The new scheme has been succesfully validated on some test 

problems for the Baer-Nunziato model of compressible multiphase flows. 
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1. Introduction 

Physical plienoinena governed by non-conservative liyperbolic systems 
arise in many scientific and technological areas, such as aerospace and au- 
tomotive industry, geophysical flows, compressible multi-phase flows, oil and 
gas extraction, among many others. The mathematical modelling of these 
phenomena is very complicated and so far no universally accepted model ex- 
ists. A common problem in many mathematical models used in the flelds 
listed above is that the governing PDE system can not be written in fully 
conservative form. A special problem that concerns compressible multi-phase 
flows is the accurate resolution of material interfaces over long times. Among 
possible strategies to achieve a good resolution of material interfaces are the 
use of Lagrangian methods [611 [HI [HI [201 [271 iH IM] , ghost-fluid and level- 
set methods [HI HSl [721 [SIl HZ] , little dissipative Riemann solvers combined 
with high order schemes [HSl [HOI [10] and, last but not least, the use of adap- 
tive mesh reflnement (AMR). Therefore, in this paper we suggest to combine 
high order accurate WENO flnite volume schemes with AMR for the solution 
of compressible multi-phase flow problems to assure an accurate resolution 
of the material interfaces. 

One particular and rather widespread mathematical model used in multi- 
phase flow applications is based on the deflagration-to-detonation transition 
(DDT) theory, from which the Baer-Nunziato equations [7] can be derived. 
Reduced models have been proposed in [551 [HHl [93] , where the original Baer- 
Nunziato system has been simplifled by carrying out an asymptotic analysis 
in the limit of stiff relaxation source terms. 

Other models in the field of multi-phase fiows are those proposed and 
analyzed by Saurel and Abgrall [THl [5l [791 [HE] for the mixture of two com- 
pressible fluids, the depth-averaged debris-flow model by Pitman and Le [75] 
as well as single and multi-layer shallow water equations [211 [731 13] ■ 

All the models above can be cast into the following general form of a 
nonlinear system of PDE in multiple space dimensions 

^ + V-F(u) + B(u)-Vu = S(u), (1) 

where u is the state vector; F(u) = [f(u), g(u), h(u)] is the flux tensor for 
the conservative part of the PDE system, with f(u), g(u) and h(u) ex- 
pressing the fluxes along the x, y and z directions, respectively; B(u) = 
[Bi(u),B2(u), B3(u)] represents the non- conservative part of the system, 



written in block- matrix notation. Finally, S(u) is the source term, which 
may in principle be stiff. When written in quasilinear form, the system ([I]) 
becomes 

^ + A(u).Vu = S(u), (2) 

where the matrix A(u) = [Ai, A2, A3] = c?F(u)/9u + B(u) contains both 
the conservative and the non-conservative contributions. 

Recent work on numerical schemes for systems of equations involving 
non-conservative terms, like Eq. (jT|, includes the family of so-called path- 
conservative schemes [211 SSI ES, d] which are based on the theory proposed 
by Dal Maso, Le Floch and Murat [65] and are a generalization of the usual 
concept of conservative methods for systems of conservation laws. Note that 
the weak formulation of the Roe method by Toumi jHI] can also be con- 
sidered as a path-conservative scheme. It has to be clearly stressed that 
path-conservative schemes have known deficiencies, which have been studied 
in detail in HES]. 

In this paper we provide the first implementation of high order path- 
conservative schemes for non-conservative systems of the type ([I]) using the 
ADER approach together with space-time Adaptive Mesh Refinement (AMR). 
In this respect, the present work can be considered an extension of the method 
proposed in [42], where an ADER-WENO AMR scheme was proposed for 
the conservative case. ADER schemes, originally developed by E.F. Toro 
and collaborators in [SH| and extensively used in the context of hyperbolic 
problems [851 ESI EH El E], are a class of numerical methods that obtain 
high order one-step schemes in time without the use of backward time levels, 
like in Adams-Bashforth type time integrators, and also without the use of 
substages, as used inside Runge-Kutta time integrators. The feature of high 
order one-step time integration will be actually the key for the construction 
of reasonably simple high order accurate AMR schemes together with time 
accurate local time stepping (LTS). For previous works including LTS see, 

e.g., [IHllMllHIlEallHniESlEaEHlllS]. 

However, the original ADER approach as proposed in [891 ESI E3] suf- 
fers from the drawback that it uses Taylor expansions in time where time 
derivatives are then replaced by spatial derivatives through a repeated use 
of the governing system of equations. This procedure, also known as the 
Cauchy-Kowalewski procedure, becomes rather cumbersome when dealing 
with complex systems of equations, and it fails completely in the presence 
of stiff source terms. A successful alternative to this strategy, that we also 



follow in this work, was proposed in [3S]. There, a different formulation of 
the ADER approach was developed, where the Taylor series expansions and 
the Cauchy-Kovalewski procedure have been replaced by a local space-time 
Galerkin method, i.e. by a weak formulation of the PDE in space-time, 
see also the references [3Sl SH El IH] where this new version of the ADER 
approach has been used. 

As for the AMR aspects of our work, we have followed a cell-by-cell refine- 
ment strategy, which is particularly convenient within our high order one-step 
finite volume approach. Various examples of AMR schemes can be found in 
literature, for an overview see the original work by Berger et al. [T71 HSl [IH 
Ml [13] and other implementations inPHUlESlEniinilZSlEIllHlESlESllS]. 

The outline of the paper is as follows: In Sect. [2]we present our numerical 
method, with some emphasis on the implementation of path-conservative 
schemes within the ADER approach. Sect. |3] is devoted to the description 
of the Adaptive Mesh Refinement infrastructure, while in Sect. |4]the Baer- 
Nunziato equations are recalled. Sect. [5] contains a set of numerical test 
problems and computational results to validate the proposed high order path- 
conservative ADER-WENO AMR schemes, and, finally. Sect. [6] reports some 
conclusions of our work and possible future extensions. 

2. Numerical method 

In order to obtain a numerical solution of the problem ([I| we use higher 
order finite volume methods in the context of the ADER framework. To 
simplify the presentation, we first describe all details of the algorithm for a 
uniform Cartesian grid. The AMR technique will be described later. We re- 
call that within the finite volume methodology, the numerical solution of the 
evolved quantities is represented at the beginning of each time-step by piece- 
wise constant cell averages. The update of these data and the computation 
of the corresponding numerical fluxes can be performed using higher order 
piecewise polynomials of degree M that have to be reconstructed, starting 
from the underlying piecewise constant cell averages, see jSU [3H | 



2.1. The Finite Volume scheme for nonconservative system,s 

The system dl]) is written in Cartesian coordinates and in three space 
dimensions as 

dvL di dg dh r. 9u -o du r, 9u 

dt dx dy dz dx dy dz 



In order to obtain the finite volume representation of pi), we discretize the 
computational domain Q in space-time control volumes defined as Xijk = 
hjk X [r,t" + At] = [x,_i,x,+ i] X [2/,_i,2/,-+i] X [Zf,_i,z^^] X [t",t" + At], 
with Axi = x^+i - Xi_i, Ayj = y-^ - y._i_, Azk = z^^ - z^.i and 
At = t"+^ — t". Each space control volume lijk defines a computational cell, 
which will be denoted so forth as Cm, identified by its mono-index m, with 
1 < m < A'^e, where Ne is the number of computational cells in the domain. 
After integration of (tsl) over a space-time control volume Xijk one obtains the 
following finite volume formulation: 
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is the spatial average of the solution in the element lijk at time t", while 
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are the average fluxes along each Cartesian direction. Furthermore we have 
defined the space-time average of the smooth part of the non-conservative 
product as 

P.,, = ^^^^ I I I I ^^""'^ ■ ^""^ ^' ^^ ^^ ^' ^^^ 

*" ^,_i Vi-i ^fc_i 

and the space-time averaged source term 

,n + l-,+ l!/,+ l-.+ l 

S.,. = i^^^^| / / J Sici,ix,y,z,t))dzdydxdt. (10) 

* 2 J 2 2 

The terms q/i in Eq. (|6|-(10) are piecewise space-time polynomials of degree 
M and represent the time-evolved reconstruction polynomials. To obtain 
the q/i first a WENO reconstruction polynomial w^ is obtained from the cell 



averages u^k at time t" (see Sect 2.2) and subsequently the time evolution 



is carried out via a local space-time DG predictor as illustrated in Sect. 2^ 
In order to integrate the non-conservative product we use the Dal Maso- 
Le Floch-Murat theory (see [65]) where the non-smooth part of the non- 
conservative term is defined as a Borel measure. In this formulation we 
therefore also need to account for the jumps of q^ at the element boundaries 
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using the path integrals 



^ B„(vl/(q„q+,.))^d., (12) 



where \l/(s) is a path joining the left and right boundary extrapolated states 
q^ and q^ in state space. The simplest option is to use a straight-line 
segment path 



^ = ^(Qh ' q^' ^) = q/. + Mt - ^h ] 



< s < 1. 



(13) 



Though simplified, this particular choice is useful in many applications since 
in the case of the shallow water equations it guarantees that the resulting 
numerical scheme is well-balanced and for the Baer-Nunziato model it has 
been shown to preserve the Abgrall condition [H |2] when used with FORCE 
and Osher-type Riemann solvers [361 HO]. This may be no longer the case for 
other systems of equations, for which eventually more sophisticated paths 



must be adopted, see e.g. [68]. With the choice of the path (13), the terms 



Dm in (12) can be computed as 



V^iq-,qt)=(^l^ B^(^(q-,q+,.))ci.^ [qt - q^) ■ (14) 



The practical computation of the integrals ( 14 ) is typically performed through 



a three-point Gauss-Legendre formula [36, 39lll0]. Finally, for the numeri- 
cal approximation of the fiuxes (l6])-(J8]) we have either adopted a local Lax- 
Friedrichs fiux (Rusanov fiux) or a simplified Osher-Solomon fiux formula 
proposed in [3211101120], 



H^H,^t) = l{Kcii 



f(q.)) 




A^inds (q^-q^:), (15) 



with 



lAil = R|A|R" 



|A| = diag(|Ai|,|A2 



|A^|) 



(16) 



and where the path \E' in Eq. (15) is the same segment path adopted in (13) 



for the computation of the jumps Vm- Again, the path integral is evaluated 
numerically using Gauss-Legendre quadrature rules. An entirely analogous 
procedure allows the computation of the numerical fiuxes g and h. Note that 



in the formulation above the numerical fiux (15) contains in its dissipative 
term both, the conservative and the non-conservative part of the system (|I|. 



2.2. WENO reconstruction 

In order to compute high order intercell fluxes and to integrate the source 
terms, it is necessary to carry out a reconstruction from the cell averages of 
the solution available at the beginning of each time step. To this extent, we 
provide the basic information about the WENO implementation that we have 
adopted, which differs from the standard one by [54j- Although a genuine 
multidimensional reconstruction is very natural within finite volume methods 
and is applicable to grids with arbitrary triangulations, we have adopted a 
cheaper dimension-by-dimension methodology [Mf ES], which is also much 
simpler to implement in the presence of adaptive mesh refinement. 

A crucial aspect to consider is the choice of the basis functions to be used 
in the reconstruction process. Two alternative options are available: a modal 
basis or a nodal basis, both of them rescaled on a reference unit interval, 
e.g. I = [0; 1], through the following changes of coordinates valid for each 
element lijk 

e = ax,^) = ^{x-x,_.) , (17) 



V = v{yJ) = ^ (y-Vj-^J , (18) 

C = a^,k) = -^(z-z,A . (19) 



Azk 



2 



The modal basis is formed by a set of M -|- 1 linearly independent poly- 
nomials, typically the Legendre polynomials, having degree from zero to M. 
The nodal basis is instead formed by a set of M + 1 linearly independent 
polynomials, {V'/}/=i'^) all of degree M, which are effectively the Lagrange 
polynomials interpolating a set of M + 1 nodal points, {xk}liX^, in such a 
way that 

tPi{xk) = 5ik l,k = l,2,...,M + l. (20) 

Numerical experiments carried out by [53] have shown that it is more efficient 
to use the nodal basis instead of the modal one, especially if the Gauss- 
Legendre nodes are used. The reconstruction is performed for each cell lijk 
on a reconstruction stencil that, for each Cartesian direction, is given by 

i+R j+R k+R. 

e=i—L e=j—L e=k—L 



where L and i?, which depend both on the order and on the specific stencil 
considered, are the spatial extension of the stencil to the left and to the right, 
respectively. Odd order schemes (even polynomial degrees M) always adopt 
three stencils, one central stencil (s = 1, L = i? = M/2), one fully left-sided 
stencil (s = 2, L = M, i? = 0) and one fully right-sided stencil (s = 3, L = 0, 
R = M). Even order schemes (odd polynomial degree M) always adopt four 
stencils, two of which are central (s = 0, L =floor(M/2) + 1, R =fioor(M/2)) 
and {s = 1, L =fioor(M/2), R =fioor(M/2) + 1), while the remaining two 
are again given by the fully left-sided and by the fully right-sided stencil, as 
defined before. The total amount of cells of each stencil is the same as that 
of the order of the scheme, namely M + 1. 

2.2.1. Reconstruction along X 

The reconstruction is first performed in the x direction, by writing the 
reconstructed polynomial in terms of the nodal basis ipi{C) 

M 

where we have used the Einstein summation convention, implying summa- 
tion over indices appearing twice. Integral conservation on all elements of 
the stencil then yields the linear equation system from which the unknown 
reconstruction coefficients w"'f can be determined: 






^e Jx 1 



Once the reconstruction has been performed for each of the stencils relative to 
the element lijk, we finally construct a data-dependent nonlinear combination 
of the polynomials obtained for each stencil, i.e. 



Ns 



wUx, t") = MO^?jk,p, with wi;.,^p = J2 ^s^7fk,p^ (24) 



s=l 



where, as specified at the beginning of Sect. |2.2[ the number of stencils is 
Ns = ^ 01 Ng = 4, for even or odd M, respectively, while the nonlinear 
weights are given by the relations 

22a ^q (^s + e) 



The oscillation indicator cr^ is 

as = Sz^wf'^w:^^ , (26) 

and it requires the computation of the oscillation indicator matrix (see |35j ) 

^Im- 2^J —^^ W^"^^ ■ ^^^^ 





^ J (9^" 9e 



Unlike the original pointwise WENO of [51] , which is 2M + 1 order accurate 
in smooth regions of the solution, our M + 1 order accurate implementation 
of WENO allows for a pragmatic choice of the coefficients A^. In particular, 
we select As = 1 for the one-sided stencils and As = 10^ for the central 
stencils. Moreover, we use e = 10~^^ and r = 8. 

2.2.2. Reconstruction along y and z 

Because the resulting reconstructed polynomial w^(a;,t") is only a poly- 
nomial in X direction, but still an average in the y and z directions, the 
reconstruction algorithm described above must be applied again along the 



remaining two directions. In practice, the steps from (22) to (24) are re 



peated and details about this procedure can be found in ^ 



2.3. The local space-time Galerkin predictor 

The objective of the local space-time Galerkin predictor method is to 
provide the time evolution, locally for each element, of the reconstructed 
polynomials w/i(x, t") obtained through the WENO reconstruction described 
before. However, unlike the original ADER approach of Titarev and Toro, 
which was based on a Taylor expansion in time and required a repeated use 
of the governing PDE system in order to substitute time derivatives with 
space derivatives, the new method relies on a weak integral formulation of 
the governing PDE in space-time using an element-local space-time Galerkin 
method. As a result, all that is required is a point-wise evaluation of fluxes 
and source terms. The result of the local space-time Galerkin predictor are 
the high order space-time polynomials q/j that are needed for the evaluation 
of the numerical fluxes, source terms and nonconservative jump terms in the 
scheme Q. In the following we briefly illustrate the method, referring to 
m ES Ellia for more details. 
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After introducing the reference time coordinate r = (t — t")/At, we write 
the system (llj) in terms of the reference coordinates r and ^ = {^,r],Q, see 
the definitions (17)-(19), to get 
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(29) 



(30) 



Ax, ^' ' Ay, ^' ^ A^fc 

To obtain a local space-time discontinuous Galerkin approach, we multi- 
ply expression (28 ) by piecewise space-time polynomials ^q(^, rj, (, r) of degree 
M, which are given by a tensor-product of the basis functions ipi adopted in 
the reconstruction procedure. Here, we use the multi-index q = {p,q,r,s). 
Integration over the space-time reference control volume then yields 
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B;^ ] d^dr]dCdT, (31) 



The discrete space-time solution of Eq. (31) will be henceforth denoted by 
q/i, and is expanded over the same space-time basis of polynomials as 

qh(e>) = ^p(e>)qp, (32) 

where qp are the unknown nodal degrees of freedom. A similar nodal repre- 
sentation is provided for the remaining terms entering Eq. (31). Integration 
by parts in time yields 
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Eq. (33) should be regarded as a nonlinear algebraic equation, to be solved 
locally for each element of the computational grid in the unknowns qp. Ad- 
ditional details about the local space-time Galerkin predictor for the specific 
case of non-conservative systems can be found for example in |36] . We em- 
phasize that the choice of a nodal basis based on Gauss-Legendre nodes allows 



a dimension-by- dimension evaluation of the terms appearing in Eq. (33). 



3. Adaptive Mesh Refinement 

A detailed illustration of the AMR implementation within our ADER- 
WENO approach has been presented in jl2]. The description we provide 
here is self-contained, but focused on the essential aspects only. Unlike the 
original patch-based block-structured approach by Berger & Oliger [T71 dSl 
[H], we refine individual Cartesian cells, which are treated as elements of 
a tree data structure, like in [57] . This choice is particularly suited to our 
element-local space-time DG predictor, which does not need any exchange of 
information through neighbor elements, and can therefore be implemented 
with no modifications even if two adjacent cells belong to two different levels 
of grid refinement. 

3.1. AMR implementation 

Each level of refinement is indicated with i, ranging from the coarsest 
level i = to the maximally refined level i = ima,x, beyond which no further 
refinement is possible. In addition, we use Ci to denote the union of all 
elements up to level i. Any cell Cm, at any level of refinement, is identified 
with a unique positive integer number m, with 1 < m < Nqciis, where A^ceiis 
is the (time-dependent) total number of cells at any given time. When a 
cell Cm at level i is refined, we refer to it as a mother or parent cell, and 

12 



the cells on the next refinement level i + 1 contained in it are called children 
cells. Moreover, the Neumann neighbors Mm of a cell Cm are the neighbor 
cells that share a face with Cm- Each cell has 2d Neumann neighbors in d 
space dimensions, except of course the case of the cells at the boundaries of 
the computational domain. On the other hand, the Voronoi neighbors Vm of 
a cell Cm are those cells which share at least one node with Cm-, and each cell 
has S'^ — 1 Voronoi neighbors in d space dimensions. 

Since the high order finite volume schemes used in this paper need infor- 
mation from neighbors to carry out the WENO reconstruction (even more 
than just the Voronoi neighbors Vm), each refined cell on a level £ + 1 must 
be surrounded by a layer of either real or virtual cells on the same level. 
The layer thickness must be greater or equal to the size of the reconstruction 
stencil. Likewise, a mother cell on level £ that is refined continues to exist as 
a virtual mother cell since it may be surrounded by non-refined cells on the 
same level which need its information for reconstruction. A schematic repre- 
sentation of this mechanism involving one level of refinement is reported in 
figure [l} There, the central cell of level £ is refined, hence it becomes virtual 
and has real children on level f + 1, while the surrounding cells are virtually 
refined in order to allow the real cells on level £ + 1 to perform the WENO 
reconstruction. 

Having introduced this terminology, it is convenient to list schematically 
the rules that are adopted in our cell-by-cell AMR implementation. 

• Any AMR scheme requires a criterion for deciding whether a given 
cell Cm needs refinement or recoarsening. We have adopted the same 
strategy described in [62] , which is based on the calculation of a second 
derivative error. A cell Cm is marked for refinement if Xm > Xret, while 
it is marked for recoarsening if Xm < Xrec, where 



Xr. 



j:,i{d^<^/dxkdxiy 



Y.,i[{\d<^/dxkU+i + \d<!>/dxk\,)/Axi + e\{dydxkdxi)m]^- 

(34) 
The summation Ylk i ^^ taken over the number of space dimension of 
the problem in order to include the cross term derivatives, whereas 
$ = $(u) is a generic indicator function of the conservative variables 
u. In most cases we have adopted Xref in the range ~ [0.2, 0.25] and 
Xrec in the range ~ [0.05, 0.15]. Finally, the parameter e acts as a filter 
which prevents refinement in regions of small ripples and is given the 
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Figure 1: Sketch of a cell- by-cell refinement of the central cell Cm on level £. The children 
of Cm are real cells on level £ +1, surrounded by virtual children of the Voronoi neighbors 
Vm of cell Cm- The virtual children are on the same level £ + 1 and are needed for 
reconstruction. 



value e = 0.01. 

Whenever a mother cell of the level 
cells, such that 



is refined, it generates x'^ children 



Axi = tAxi+i Aye = rAy^+i Azi = xAzi+i. 



(35) 



In addition, and as commented below, the time steps are also chosen 
locally on each level so that 



Atf = xAt 



i+i- 



(36) 



Due to the high order WENO reconstruction, the refinement factor t 
must satisfy t > M. 

At any level of refinement, each cell Cm has one among three possible 
status fiags, which we denote by a. The first possibility is that Cm is an 
active cell (a = 0), in which case it is updated through the finite volume 
scheme described in the previous Section [2j The second possibility is 
that Cm is a virtual child cell (a = 1) and is updated by projection 
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of the mother's high order space-time polynomial. In practice, the 
virtual children receive their values from the active mother via standard 
L2 projection. For this purpose, the space-time polynomials qh can 
be conveniently evaluated at any time. This operation is needed for 
performing the reconstruction on the finer grid level at intermediate 
times. The projection operator for a cell Cm on level i is simply given 
by evaluating the space-time polynomial q/i of its mother at any given 
time t^ as follows: 

Finally, Cm can be a virtual mother cell {a = —1), updated by recur- 
sively averaging over all children from higher refinement levels. Namely, 
the virtual mother cell obtains its cell average by averaging recursively 
over the cell averages of all its children, namely including the possible 
children of their children. If we denote the set of children of a cell Cm 
by Bm, then the averaging operator is given by 

• Only real cells [a = 0) can be refined. Therefore, if a virtual cell needs 
to be refined it must be first activated. 

• The levels of refinement of two cells that are Voronoi neighbors of each 
other can only differ by at most unity. Moreover, every cell has Voronoi 
neighbors, which can be either active or virtual, at the same level of 
refinement. 

3.2. AMR local time stepping 

As anticipated above, each of the refinement levels is advanced in time 
with its own local time-step, i.e. At^ = rAt^+i. The use of time steps that 
are integer multiples of each other among subsequent levels is very conve- 
nient, and indeed very natural within AMR. However, alternative local time 
stepping schemes are also possible (see [Ml EB E3]), where a different local 
time step is allowed for each element. After denoting by t" and t?"*"^ the 
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Figure 2: Example of a local time stepping algorithm involving two levels of refinement 
with r = 2. Legend: • (regular active cell), O (virtual refined cell), D (virtual coarse cell). 
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current and future times of the level £, the first level to be updated is the 
one with the largest value of £ satisfying the update criteriorPl [38] 

t1^'<ttl 0<£<Cax. (39) 

In practice, starting from the common initial time t = 0, the finest level of 
refinement £max is evolved first and performs a number of r sub-timesteps 
before the next coarser level £max ~ 1 performs its first time update. This 
procedure is then applied recursively and it implies a total amount of r^ sub- 
timesteps on each level to be performed in order to reach the time tg^^ of 
the coarsest level. As example of the application of the local time stepping 
strategy to a one dimensional case involving two levels of refinement with 
r = 2 is reported in Fig|2] 

Thanks to the use of the local space-time predictor, which computes the 
predictor solution (\h for each element valid from time t" to time t""*"^, the 
computation of numerical fluxes between two adjacent cells on different lev- 
els of refinement is rather straightforward. Further details about the actual 
implementation of the local time stepping procedure and of the AMR par- 
allelization through the standard Message Passing Interface (MPI) can be 
found in |l2| . 



4. The Baer-Nunziato equations 

A Baer-Nunziato type model for compressible two-phase flow without 
relaxation terms is given by the following system of equations, see [3 EHl El 
691: 



9. 

at 
at 

at 

a_ 
at 



(0iPi) + V ■ (0ipivi) = 0, 
(0iPiVi) + V ■ (0ipivivi) + V0ipi = P/V01, 
{(pipiEi) + V ■ {{(pipiEi + (pipi) vi) = -pjdtcpi, 
{(P2P2) + V ■ (02P2V2) = 0, 



^ (02P2V2) + V ■ (02P2V2V2) + V02P2 = Pi'^h, 

at 



ai^i 



'2P2E2) + V ■ ((02P2^2 + (P2P2) V2) = PidtCpl, 
+ VjV0i = 0. 



(40) 



We define t"i[ := ig , so that also the scheduUng of level i ~ is taken into account. 
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The system is closed by the stiffened equation of state (EOS): 

Pkllk - 1) 

Here, 0^ denotes the volume fraction of phase fc, p^ is the density, v^ is the 
velocity vector, Ek = ek + |v^ and e^ are the phase specific total and internal 
energies, respectively. In the literature, one of the phases is often also called 
the solid phase and the other one the gas phase. Defining arbitrarily the 
first phase as the solid phase in the rest of the paper we will therefore use 
the subscripts 1 and s as well as 2 and g as synonyms. For the interface 
velocity u/ and pressure pi we choose v/ = vi and pi = P2, according to [7], 
although other choices are also possible (see e.g. the paper by Saurel and 
Abgrall [78j). The state vector u is 

U = (0lPl, 0lPlVi, (plPlEi, 02P2, 02P2V2, 02P2-E2, 0l) . (42) 

The system (40 ) can be cast in the form prescribed by ([I]) by collecting all the 
non-conservative terms in the matrix B(u), while keeping the conservative 
part of the system expressed through F(u). An exhaustive treatment of 
the mathematical properties of the Baer-Nunziato equations can be found in 
[SI [29|, Eni |87] , where also exact and approximate solutions to the Riemann 
problem are given. 

5. Test Problems 

In all test problems shown below the indicator function for the refinement 
and recoarsening criterion has been chosen as 



with some reference densities pi o and p2,o, respectively. 

5.1. Smooth vortex problem 

The first test that we have considered is given by a stationary and axisym- 
metric solution of the Baer-Nunziato equations and has first been reported 
in [3S]- The resulting configuration describes a vortex-type solution with no 
radial motion. Because of these assumptions, the continuity and the energy 



equations are automatically satisfied. After choosing a simple dependence 
on the radius r of the pressure and of the volume fraction of the solid phase, 
namely 

Pt = mo(l-je(l-''V4)^, (t = l,2), (44) 

o Z\J 111 

the momentum equations can be easily solved, to provide the velocity field 
of the vortex as 



2siD 

^^/2 

2p2S2 



IrD 



Pio (^4V2^Fi + 6Hi - 12Gsl + 3Hislj + 3p2o4 (4G - ^2)] (46) 
\/P2P20F2 , (47) 



where 



Hk = e 24 , Fk = e 4 , (fc = l,2), (48) 

and 

G = e-"'/^ I^ = Pi (2^2^ + 3G] . (49) 

In order to make the test problem unsteady, the vortex is then boosted along 
the diagonal of the computational domain through a Galilean transformation 
of the velocity, with components u = v. In our tests we have chosen the 
following parameters 

3 3 7 

Pi = 1, P2 = 2, pio = 1, P2o=2' *i = 2' *^ " 5' "" = ^ = 2, 

(50) 
while the computational domain is i7 = [—10; 10] x [—10; 10] with four pe- 
riodic boundary conditions, in such a way that the exact solution of the 
problem is given by the initial condition after T = 10. In Table [T] we have 
reported the results of the convergence tests, where we have used the third 
and fourth order version of the method. Here pi^o = P2,o = 1 have been cho- 
sen. The convergence rates have been computed with respect to an initially 
uniform mesh, as proposed by Berger and Oliger in [T7] . 
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Table 1: Numerical convergence results for the vortex test using the one-step ADER- 
WENO finite volume scheme. The error norms refer to the variable pi at the final time 
T — 10, and have been computed with -^max — 2. The asterisk * refers to a uniform grid. 



Ng^Ng 


eL2 


0{L2) 


Ng^Ng 


<^L2 


0{L2) 






C3 






C4 


15x15* 


4.9627E-01 
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2.3166E-02 
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Figure 3: AMR grid of the isentropic- vortex test at the initial time (left panel), at time 
t — 2.0 (central panel) and at the final time ^ = 10 (right panel). Two levels of refinement 
have been adopted (^max — 2), starting from a uniform 45 x 45 grid. 



5.2. Riemann problems 

The high order space-time adaptive ADER-WENO methods proposed in 
this paper are particularly well suited for the accurate resolution of mate- 
rial interfaces in multi-fluid and multi-phase flow problems. This claim is 
validated in the following by solving a set of ID shock tube problems on 
space-time adaptive Cartesian meshes in 2D. The exact solutions for the 
Riemann problems solved here have been provided in [6l |80l [29j. From these 
papers a set of six Riemann problems has been chosen, see Table [2j The 
same set of problems has already been solved with high order unstructured 
one-step WENO finite volume schemes using a centered path-conservative 
FORCE method in [36]. A subset of these Riemann problems has been 
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solved again in |3D] using the more accurate path-conservative extension of 
the Osher method, which is also used in this paper. 

The two-dimensional computational domain is fi = [—0.5; 0.5] x [0;1], 
which is discretized at the level grid with only 100 x 10 cells. A maxi- 
mum number of two refinement levels (i'max = 2) is chosen, together with a 
refinement factor of r = 4. The discontinuity is located initially at a; = 
and the final simulation times are given in Table [2j For all test cases a third 
order ADER-WENO scheme is used with reconstruction in characteristic 
variables. In x-direction transmissive boundary conditions are imposed and 
periodic boundaries are applied in ?/-direction. The parameters for the refine- 
ment criterion are pi.o = P2,o = 1 apart for RP2 and RP4, where pi^o = 1000 
and p2,o = 1- 

The results for the AMR computations are shown in Figs. |4]-[9j A sketch 
of the AMR grid at the final simulation time is depicted on the top left of 
each figure together with a cut through the reconstructed numerical solution 
Wfi on 200 equidistant points along the x-axis in the remaining sub-figures. 
For RP4 the same quantities as in [80] are shown. For all six problems we 
obtain an excellent agreement between the high order AMR computations 
and the exact reference solutions provided in |6|, [801 EH]- The agreement 
is much better than the one obtained in the previous publications of the 
authors [321 SD]- The solid contact is resolved perfectly well in all cases. 
Furthermore, no spurious post shock oscillations as reported in [3S] for RP5 
are visible in the present high order ADER-WENO simulations with AMR. 
Our results clearly confirm that the combination of adaptive mesh refinement 
(AMR) with high order WENO finite volume schemes with the little diffusive 
Osher-type Riemann solver [SHI HD] is very well suited for the simulation of 
compressible multi-phase flows, as claimed at the beginning of this section. 

5.3. Explosion problems in multiple space dimensions 

An explosion problem for the Baer-Nunziato equations can be solved both 
in two and in three space dimensions, after setting-up the following initial 

conditions 

, . ( Ui a r < R, , , 

u(x, 0) = < .r „ (51) 



where x is the vector of spatial coordinates, r = v x^, and R is the radius of 
the initial discontinuity, which we have set equal to 0.5 and to 0.4 for the two- 
dimensional and for the three-dimensional explosion tests, respectively. The 
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Figure 4: Results for Riemann problem RPl. 
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Figure 5: Results for Riemann problem RP2. 
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Figure 6: Results for Riemann problem. RP3. 
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Figure 7: Results for Riemann problem RP4. 



25 



''^^^'' 




■ '■ '■ '■■■■> .. .i... .1. r 11 



■n.s -cij -03 -oj -oli fl 0.1 en 03 a4 05 



1 1 - 

uu - 
1XK - 

IjM - 

Iffi - 



««& - 

*W - 
afl2 - 



<>* 



_!_ J 1 1 I J 



S -0 J -0_3 -0^ OlI C 1 a ? 05 O & 5 



^.5. -flj 



_l I I I 



•dj -03 -o_9 dLi o 0.1 aj 03 a4 ns 



^^ 



_J -i 1 1 I J 



as 


F- 


- 


ji-nr 




-02 

3- -OA 

'OS 


ADf A-WENO 01 (i 


■ 




/ 


** 









B -0 J -0.3 -OJ OlI C 1 ? 5 at. 5 



'□B -Ojt -03 -OJ -OlI O 0.1 OJ 03 Ol4 05 



1 IS 






EudulijIiiH 


jUH} 








1 1 




ACEB-WEHO Ca 




1.0S 


- 


















c£ 1 














, 










\ 












Ol9 


■ 




1 1 1 












B 


-QJ 


-0.3 -OjJ OlI 


c 


01 


03 


05 


Od 


0.5 



IJ 


|- 






nil) 


1 


. 


Jtf)M-WeHOOJ(i 


OJ 


■ 




/ 






j— 




(M. 


. ' 








_il 




*? 







iris -OjI -03 -0.3 OlI O 0.1 OJ 03 ft* 05 



Figure 8: Results for Riemann problem RP5. 
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Figure 9: Results for Riemann problem RP6. 
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inner and outer states Uj and Uo for the two test cases that we have consid- 
ered are reported in detail in Table |3| Analogously to the compressible Euler 
equations, a reference solution can be obtained after solving an equivalent 
one-dimensional problem in cylindrical coordinates for the two-dimensional 
explosion, and in spherical coordinates for the three-dimensional one. This 
essentially implies that the equivalent one-dimensional problem contains ad- 
ditional algebraic source terms on the right hand side of the equations, which 
account for the use of curvilinear coordinates, see [881 ES] . The ID reference 
solution has been computed using a classical second order TVD method with 
the Osher Riemann solver [311 SO], using 5000 grid cells. For EPa the pa- 
rameters for the refinement criterion are pi o = P2,o = 1 ^'Hd for EPb we use 
pifi = 1000 and p2,o = 1- 



Table 2: Initial states left (L) and right (R) for the Riemann problems solved in 2D. Values 
for 7i, TTi and the final time ie are also given. 





Ps 


Us Ps 


P9 


Ug 


Pg 


0s 


te 


RPl [29]: 




7s = 1.4, TTs 


= 0, 


7. = 1-4, 


TTg = 






L 
R 


1.0 
2.0 


0.0 1.0 
0.0 2.0 


0.5 
1.5 


0.0 
0.0 


1.0 
2.0 


0.4 
0.8 


0.10 


RP2 [29]: 




7s = 3.0, TT, = 


= 100, 


7, = 1.4, 


7r, = 






L 
R 


800.0 
1000.0 


0.0 500.0 
0.0 600.0 


1.5 
1.0 


0.0 
0.0 


2.0 
1.0 


0.4 
0.3 


0.10 


RP3 |29i: 




7s = 1.4, TTs 


= 0, 


7. = 1.4, 


7r, = 






L 
R 


1.0 
1.0 


0.9 2.5 
0.0 1.0 


1.0 

1.2 


0.0 
1.0 


1.0 
2.0 


0.9 
0.2 


0.10 


RP4|80]: 




7s = 3.0, TTs = 


3400, 


Ig = 1.35 


^9 = 







L 
R 


1900.0 
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0.0 10.0 
0.0 1000.0 


2.0 
1.0 


0.0 
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3.0 
1.0 


0.2 
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0.15 


RP5 m-- 




7s = 1.4, TTs 


= 0, 


7. = 1.4, 


7r, = 






L 
R 


1.0 
1.0 


0.0 1.0 
0.0 1.0 


0.2 
1.0 


0.0 
0.0 


0.3 
1.0 


0.8 
0.3 


0.20 


RP6 [5]: 




7s = 1-4, TTs 


= 0, 


7. = 1.4, 


7r, = 






L 
R 


0.2068 
2.2263 


1.4166 0.0416 
0.9366 6.0 


0.580C 
0.489C 


) 1.5833 
) -0.70138 


1.375 
0.986 


0.1 
0.2 


0.10 
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5.3.1. 2D computations 

We have first evolved the two models EPa and EPb in two spatial di- 
mensions, by adopting a fourth order ADER-WENO-AMR scheme with two 
levels of refinement over a computational domain given by [—1; 1] x [1; 1]. 
The level zero grid is composed by 50 x 50 cells, which are immediately re- 
fined according to the refinement criterion applied to the initial conditions. 
A representative example of the grid at time t = and at the final time 
t = 0.2 is shown in Fig. 10 for the model EPb. The final grid (right panel) 



is composed by 50020 cells. Fig. 11 and Fig. 12 report the results of the 
computation for the two models EPa and EPb by comparing them with the 
reference solution. 



5.3.2. 3D computations 

We have then evolved the same models EPa and EPb in three spatial 
dimensions, by adopting a third order ADER-WENO-AMR scheme with two 
levels of refinement over a computational domain given by [—1; 1] x [1; 1] x 
[—1; 1]. The level zero grid contains 34 x 34 x 34 cells. The final grid at 
time t = 0.15, shown as a representative example in Fig. 13 for the model 



EPb, is composed by 3, 833, 016 cells. Fig. [14] and Fig. [T5| on the other hand, 
report the solution at time t = 0.15, compared to the reference one. All the 
relevant features and waves of the solution are successfully resolved by the 
scheme, which remains essentially non-oscillatory and performs the correct 
grid refinement where this is needed. 



Table 3: Inner and outer initial states for the two multidimensional explosion test prob- 
lems. 
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Figure 10: AMR grid structure of the explosion test EPb in two space dimensions, at 
time t = (left panel) and at the final time t = 0.2 (right panel). Two levels of refinement 
have been adopted (£max = 2). 



5.4- Shock-bubble interaction 

In this section a strongly simplified shock-bubble interaction-type prob- 
lem is considered. The initial condition is given by a Riemann problem that 
leads to a strong isolated shock wave travelling with a shock speed of s = 100. 
The pressure jumps over six orders of magnitude across the shock. The initial 
discontinuity is located at x = and the left and right initial states, which 
are connected by a Hugoniot curve, are summarized in Table |4j A bubble of 
radius R = 0.25 is located at a; = 0.5, y = 0. The state inside the bubble is 
also given in Table |4| The parameters for the equation of state of each phase 
are 71 = 3.0, tti = 100, 72 = 1.4 and tt2 = 0. The problem is solved until 
time t = 0.0025 in a computational domain Q = [—0.5; 3.0] x [—0.75; 0.75] 
that is discretized on the level i = grid using 70 x 30 grid zones. Two 
levels of refinement are used (^max = 2) with a refinement factor of r = 4, 
which corresponds to a uniform fine grid resolution of 1120 x 480. Periodic 
boundary conditions are employed in y direction and Dirichlet boundary con- 
ditions corresponding to the left and right state are imposed in x direction. 
We choose pi = 1000 and p2,o = 1 ioi the indicator function $. The evolu- 
tion of the liquid density pi and the liquid volume fraction (pi are depicted 



for various times in Fig. 16 One can observe how the bubble is compressed 
and accelerated by the incident shock wave and how Richtmyer-Meshkov- 
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Figure 11: Results of the 2D explosion problem EPa at time t = 0.2. A cut of various 
quantities along the x-axis is reported, both for the solid (left panels) and for the gas 
phase (right panels). The ID reference solution is also shown for comparison. Two levels 
of refinement have been adopted (^max = 2). 
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Figure 12: Results of the 2D explosion problem EPb at time t = 0.2. A cut of various 
quantities along the x-axis is reported, both for the solid (left panels) and for the gas 
phase (right panels). The ID reference solution is also shown for comparison. Two levels 
of refinement have been adopted (^max = 2). 

32 




Figure 13: AMR grid structure of the explosion test EPb in three space dimensions. Two 
levels of refinement have been adopted (£max = 2), starting from a uniform 34 x 34 x 34 
grid. 



type flow instabilities occur at the bubble border at later times. One further 
notes the shock wave reflected by the bubble. A zoom into the AMR grid 
at the final time t = 0.0025 is depicted in Fig. 17 The results of this sec- 
tion are only considered as qualitative, in order to test the robustness of our 
high order AMR method in the presence of strong shock waves in liquids. 
Here, the high order ADER-WENO AMR scheme has been applied to the full 
seven equation Baer-Nunziato model, but without taking into account any 
interphase drag and pressure relaxation, hence no comparison with experi- 
mental data can be made for our results. A very detailed quantitative study 
of shock-bubble interactions with comparison against experimental data has 
been carried out in [77] using second order accurate high resolution shock 
capturing schemes together with AMR. 



Table 4: Left (L), right (R) and bubble (B) state for the shock-bubble interaction problem. 
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Figure 14: Results of the 3D explosion problem EPa at time t ~ 0.15. A eut of various 
quantities along the x-axis is reported, both for the solid (left panels) and for the gas 
phase (right panels). The ID reference solution is also shown for eomparison. Two levels 
of refinement have been adopted (^max = 2), starting from a uniform 34 x 34 x 34 grid. 

34 



H«>*r«Tic« 



1 


: 




a 


UlER MEMO 03 (UHj 






^ 


- RAlAmKA 


a 


: 


i 






t.e 


- 


1 






M 


- 


1 






\3 


^ 




/^ 


11 






, 





OS 


04 


m C.6 



r 




AKRWENaO 


slMlft) 


J 


[ 


1 


1 


■ ■ 


0! « 


4 it, o.s 







»6 C.6 



g M>En WENO 03 |un) 




kDEfl WEKKI 03 |Un| 




Figure 15: Results of the 3D explosion problem EPb at time t = 0.15. A cut of various 
quantities along the x-axis is reported, both for the solid (left panels) and for the gas 
phase (right panels). The ID reference solution is also shown for comparison. Two levels 
of refinement have been adopted (^max = 2), starting from a uniform 34 x 34 x 34 grid. 
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Figure 16: Liquid volume fraction (left) and liquid density (right) of the shock-bubble 
interaction problem at times t = 0.005, t = 0.010, t = 0.015, t = 0.020 and t = 0.025 
(from top to bottom) using a third order ADER-WENO scheme with AMR. 



36 



0.4 p. 

0.3 I 

0.2 

01 

-0.1 
-0.2 
-0.3 i 



-GA 



4-f ^Wi ■■ ■ ■ i¥---^"----------- 


^1 y^MOTL^ 


, — -^&iJ^^^^B| ijy — 


y ^^H^^^KCT^fcin b™ 


, iBIBB" |Hi 


I^hI 


r 1— WSt 


1 Pll|iP|B 


^---T-||ii||Wyf^ I--- 


1^ ikilBPlhiTliiiilJ 1 1 



1.2 1.3 1.4 1.5 1.6 1.7 1,8 1.9 2 2.1 



Figure 17: Zoom into the AMR grid at time t ~ 0.025 for the shock-bubble interaction 
problem. 
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5. 5. An application to complex free surface flows 



A reduced version of the Baer-Nunziato model (40) can also be used for 



the simulation of complex non-hydrostatic free surface flows, as proposed in 



The reduced three-equation version of the PDE system (40 ) has been 
obtained in 



after introducing the following simplifying assumptions: 
i) all pressures are relative with respect to the atmospheric reference pressure, 
which is assumed to be constant and zero everywhere and for all times; ii) 
the influence of the gas phase onto the liquid can be neglected, hence the 



evolution equations for the gas phase can be dropped in (40 ); iii) the pressure 
of the liquid phase is computed by the Tait equation of state [12], which is 
a barotropic EOS and therefore also allows the energy equation of the liquid 
to be removed from the governing PDE system. The Tait EOS reads 



Pi 



(52) 



where pi,po are the liquid density and the reference liquid density at stan- 
dard conditions, respectively. The compressibility of the fluid is governed by 
the constants fco aiid 7. To avoid low Mach number problems, an artificial 
Mach number of the order M G [0.1; 0.3] is chosen for typical environmental 
free surface problems. If real compressibility effects play a role, such as in 
industrial high pressure applications, the proper values ko = 3.2 ■ 10* Pa, 
Pq = 1000 kg/m? and 7 = 7 for real water can be chosen, which lead to the 
real sound speed in water of approximately c = 1500 m/s. 



Introducing the above-mentioned simplifications into system (40) yields 
the following reduced three-equation model |3T1 132] : 



I m + V ■ (#v) 



d_ 

at 



d_ 

at 



(0pv) 



f V 
V0: 



0, 

(pvv +pl)) 



4>ps, 



(53) 



0. 



The subscript 1 has been dropped to ease notation. The mass and momen- 
tum equations are fully conservative in the system above since the interface 
pressure pi = P2 = 0, while the advection equation for the volume fraction 
remains non-conservative. In (53) the state vector is Q = (0p, 0pv, 0) and 
g is the gravity vector acting along the vertical direction y, i.e. g = (0, —g) 
with g = 9.81 m/ s^. Further details and thorough validations of this model 
can be found in [3ll |32] . 
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The last test problem of this paper where we apply the high order ADER- 
WENO scheme with adaptive mesh refinement consists of a dambreak with 
successive wave impact against a vertical wall. In this test problem the 
free surface is strongly deformed and the impact against the wall also leads 
to wave breaking. The setup of this test case is taken from |2H1 SH]- The 
rectangular computational domain is defined as fi = [0; 3.2] x [0; 2.0] with 
refiective wall boundaries on all sides of the domain apart from the top, 
which is assumed to be transmissive. The initial liquid domain is Qi = 
[0; 1.2] X [0; 0.6], where = 0.999 is set. The initial velocity is zero and the 
initial density distribution is chosen such as to obtain a hydrostatic pressure 
distribution inside the liquid, see [21] for details. Outside the liquid domain, 
we set the pressure to zero and a low but finite value for the liquid volume 
fraction is chosen (0 = 10~^). We use g = 9.81 and fco = 2.0 ■ 10^. The 
level zero grid contains only 60 x 40 cells. We use two levels of refinement 
^max = 2 and a refinement factor of r = 4. This corresponds to an equivalent 
resolution on a uniform fine mesh of 960 x 640 elements. A third order ADER- 
WENO method is used together with the little dissipative Osher scheme. 
The contours of the liquid volume fraction are depicted for three output 



times in Fig. 18 The results of the present high order AMR method are 



compared with the results on a uniform fine grid. According to Fig. 18 the 



two solutions match very well. The AMR grid at the final time t = 1.5 is 



depicted in Fig. 19 and contains 47710 cells, while the uniform fine grid 
contains 614400 cells. This means that for the present test problem the use 
of an AMR technique allows to obtain almost the same results as with a 
uniform fine grid, but with 12.88 times less elements. Concerning the AMR 
overhead we also report the average CPU time per element update (EU) in 
Table |5} which is the total wallclock time divided by the total number of 
element updates, normalized with respect to the uniform grid. This figure 
indicates the overhead introduced by the AMR machinery. For the present 
test problem, using a third order ADER-WENO AMR scheme it is only 
19 %, which agrees with the data for the AMR overhead published in |12] . 
One can note how the AMR algorithm refines the grid in the vicinity of the 
free surface, which is resolved in a very sharp manner. The fiow evolution 
computed with the present method is also in good agreement with the 3D 
SPH computations performed in |16] and with the results presented in [311 



In Fig. 20 we compare the experimental data for the pressure at the 
wall [HI] with the numerical results obtained with the method proposed in 
this paper. The amplitude and the arrival time of the first pressure peak due 
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to the wave impact at the wall is captured correctly. 

6. Conclusions 

In the present paper the first better than second order accurate path- 
conservative one-step WENO finite volume scheme with adaptive mesh re- 
finement (AMR) has been presented for the solution of non-conservative hy- 
perbolic systems. The method has been applied in particular to the Baer- 
Nunziato model of compressible multiphase fiows. It has been shown via a 
numerical convergence study that the proposed numerical approach reaches 
the designed high order of accuracy in space and time. Furthermore, the 
accuracy and robustness of the method has been validated on a large set of 
test problems, for which either exact or other reference solutions were avail- 
able. It has been clearly shown that the use of AMR can lead to a significant 
reduction in grid elements and CPU time compared to simulations performed 
on fine uniform grids. Compared to our previous publications [36| HO] , where 
high order schemes have been employed without AMR, the material inter- 
faces are much better resolved using the present high order AMR technique. 
This becomes particularly evident in the results obtained for the shock tube 



problems solved in Section 5.2 



Future applications of the present high order one-step AMR methodology 
may concern the simulation of bubbles with phase transition [30] and chem- 
ically reacting compressible multiphase fiows. For the latter it has already 
been shown before in literature that AMR techniques can be very useful 
to resolve all the length scales involved in such multi-scale problems |52] . 



Table 5: Memory and CPU time comparison of the third order ADER-WENO AMR 
method and ADER-WENO on a uniform fine grid for the dambreak and wave impact 
problem of Section [5. 5| Memory consumption is measured in number of real elements and 
CPU time is normalized with respect to the wallclock time for the fine uniform mesh. The 
normalized average time per element update (EU) is given in the last row to quantify the 
overhead of the AMR framework. 

AMR Uniform ratio 



Number of real cells 

Total CPU time 
Average time per EU 



47710 


614400 


12.877 


0.1134 


1.0 


8.816 


1.1923 


1.0 


0.839 
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Figure 18: Dambreak and wave-impact problem at times t = 0.5, t = 1.0 and t = 1.5 using 
a third order ADER-WENO scheme with AMR (left) and with a uniform fine grid (right). 
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Figure 19: AMR grid for the dambreak and wave impact problem at time t — 1.5. 
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Figure 20: Pressure evolution at the wall. Comparison of experimental data with the 
computational results. 
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Other relevant mathematical models that could be treated in the framework 
outlined in this paper include the nonconservative debris flow model of Pit- 
man and Le [75] as well as single and multi-layer shallow water equations 
pH [73] |25| |3]. In particular for multi-scale tsunami wave simulations the 
combination of a high order method that allows a coarse grid to discretize 
the wave propagation in the ocean with a flne grid that allows to follow and 
resolve all the relevant details on the coastline may be useful. 
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